Counter-ion density profile around charged cylinders: the strong-coupling needle limit 
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Charged rod-like polymers are not able to bind all their neutralizing counter-ions: a fraction of 
them evaporates while the others are said to be condensed. We study here counter-ion conden- 
sation and its ramifications, both numerically by means of Monte Carlo simulations employing a 
previously introduced powerful logarithmic sampling of radial coordinates, and analytically, with 
special emphasis on the strong-coupling regime. We focus on the thin rod, or needle limit, that 
is naturally reached under strong coulombic couplings, where the typical inter-particle spacing a 
along the rod is much larger than its radius R. This regime is complementary and opposite to the 
simpler thick rod case where a' <C R. We show that due account of counter-ion evaporation, a 
universal phenomenon in the sense that it occurs in the same clothing for both weakly and strongly 
coupled systems, allows to obtain excellent agreement between the numerical simulations and the 
strong-coupling calculations. 
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I. INTRODUCTION 



C*~) "Some linear biopolymers are intrinsically stiff objects, fur- 

ther rigidified by electric charges along their backbone. This 

cipthe case of double stranded DNA, tubulin, actin, and some 
Viruses. These macromolecules may be envisioned, to first 
Approximation, as straight charged cylinders, attracting neu- 
tralizing counter-ions though a logarithmic potential. It was 
| 'fi rst realized by Onsager that this functional form is similar 
C^lo that of the confinement entropy, so that not all counter- 
"£|*pns are bound to the polymer a finite fraction only 
Remains confined in the limit of infinite system size, with no 
External boundary; above a certain critical temperature, all 
^jounter-ions evaporate. 

1 . The phenomenon of counter- ion condensation/evaporation 
jjjas been central in a consequential number of studies since 
ghc 1970s, for it is cardinal in a wealth of static and dy- 
C aa mic properties of charged polymers [IHsj], see e.g. @- 
~m| for more exhaustive references. It is governed by the 
s rj-called Manning parameter £, which is the dimensionless 
line charge of the rod: £ = qls\ defined from the the va- 
f~iency q of counter-ions, the linear charge Ae, and the Bjer- 
CHjum length Ib = e 2 / (AnekT) where e is the dielectric con- 
■^tlant of the medium, e is the elementary charge and kT de- 
^^otes thermal energy. Interestingly, the mean-field (Poisson- 
^vpoltzmann 1121) scenario of a complete evaporation of ions 
i-Hbr £ < 1 in an infinite system, and of a partial conden- 
C^ation for £ > 1 also holds beyond mean- field @,|, when 
' # ^pulombic couplings are important. Such couplings are con- 
veniently quantified by the parameter 2 = q 2 lsS,/R where 
"^Jjj? is the cylinder radius. It is defined to match its planar 
rjjounterpart 3 = 2nq 3 l B a [l(| EH, where ae is the surface 
j^harge of the colloid considered, with here a = \/(2nR). 
While £ and 3 both measure the inverse temperature, their 
scaling with 1/T differ, and their roles in the forthcoming 
analysis are somewhat asymmetric. In essence, 3 governs 
the distance to mean-field, in the sense that the Poisson- 
Boltzmann theory becomes exact for 3 — > 0, and remains 
accurate for small 3. In the present work, we will focus 
on the opposite limit of large 3 values, which defines the 
strong-coupling regime for which mean-field is invalid. One 
may naively think that the Manning parameter should also 
be large for a strong coupling approach to hold, but we will 
see that the situation is more subtle, and some quantities 
can be obtained for arbitrary values of £ provided 3 is large. 
It is important though to clearly discriminate the thin and 
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• Counter-ion with charge -qe 

Figure 1: The Wigner crystal formed for S — > oo at weakly 
curved quasi- planar cylinder (thick cylinder case, a <C R). The 
shaded region is the cell area per counter-ion at the surface. Up 
to a numerical prefactor, we have a oc \fqfo. At finite but large 
S, such an idealized configuration is met for £ ^> H 1//2 . 



thick cylinder cases, because they involve different mecha- 
nisms. To sec this, we start with a cylinder with large ra- 
dius, where "large" means that R significantly exceeds the 
typical distance a between ions when they are close to their 
ground state. The situation is locally that depicted in Fig. 
[1] and is mostl y g overned by the planar geometry physics 
studied in [l~3l- 16|. up to some curvature corrections that 
have not been studied so far. From the electro-neutrality 
requirement aa 2 oc q : we obtain that £ 3> 3 1 / 2 for thick 
rods. For the sake of the discussion, wc restrict here on 
ground state considerations, and will address thermal effects 
in detail later on. Gradually decreasing R at fixed a (with 
fixed 3 and decreasing £), one encounters typical configu- 
rations such as that depicted in Fig. [5] where a and R are 
comparable, before reaching the thin or needle case sketched 
in Fig. [31 where £ <IC \/S (a R), and that will be the 
center of our interest. The reason is that a common and ex- 
perimentally relevant way to raise Coulombic couplings in a 
soft matter system is to increase the counter-ion valency q. 
Given that 3 x / 2 /£ cx y/q/(XR), we see that this ultimately 
leads to £ <C It may be noted that the ground state, 

reached e.g. at T — > 0, corresponds to both diverging 5 and 
£ parameters, but the ratio S 1 / 2 /^ is a geometric quantity, 
independent of the temperature. In the needle limit, the 
relevant length scale to measure counter-ion distances is no 
longer a, defined as a two-dimensional quantity, but its one 



dimensional counterpart a' (see Fig. [3]). We can obtain an 
order of magnitude assuming all ions are condensed onto the 
rod, which leads tcQ a 1 ~ q/X. This is a lower bound, since 
the phenomenon of counter-ion evaporation invariably leads 
to a lower linear charge than Ae. 




counter-ions 




Figure 4: The 3D cylindrical cell model. The rod is assumed 
positively charged, and the counter-ions are therefore negative, 
with charge — qe. 



Figure 2: Artist's view of couter-ions at the charged cylinder 
for large S, with a and R of like magnitude. Here, £ is 
comparable to 
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Figure 3: Schematic representation of the ground state in the 
needle limit, that is for £ <C vS. The distance between charges 
is denoted a' and we have tt' > a > iJ. 



As alluded to above, we consider an infinitely long charged 
cylinder of radius R, within the primitive cell model (see 
Fig. [4]): point counter-ions with charge — qe are confined 
in a coaxial larger cylinder of radius D. The ions may 
equally have a small hard-core, which is immaterial here. 
We thus deal with a salt-free system. The solvent is ac- 
counted for though its uniform dielectric constant e. No 
dielectric discontinuity is considered here between the sol- 
vent and the charged cylinder. Charged species interact with 
three-dimensional Coulomb potential, varying as inverse dis- 
tance for point particles, with an additional hard-core term 
that prevents the ions from entering the charged cylinder. 

The outline of the paper is as follows. We first remind 
some relevant and known results pertaining to mean-field 
in section |nl before presenting in section lim our strong cou- 
pling (SC) analysis. For a given value of Manning parameter 
£, we work out in a first step the leading order behavior of 
the density profile when the coupling parameter 3 — > oo. 



1 It is customary to introduce the Gouy length fj, = (27tZb<t<;) — 1 . We 
then have S = Isq 2 /^ fl7| . and likewise £ = R//J, [fj. Since the 
interesting regime is for £ > 1, fj, is always the smallest length scale 
in the problem. In the needle limit where S 1 / 2 2> f, we have in 
general |i<RC«<«'. The relation between a' and q 2 ls is given 
by £, with roughly q^ls/o, 1 ~ § (taking due account of counter-ion 
condensation, which affects a', we get q 2 lg/a' = £ — 1). 



The present problem in the needle limit is a case where in 
principle and for the leading order only, the virial strong- 
coupling approach put forward by Netz and collaborators 
@, i, Hfl, HI [H [U OH should coincide with its Wigner-SC 
counterpart worked out in [H|,[l6|. However, due to a differ- 
ent treatment of counter-ion evaporation, our leading SC ex- 
pressions (SC-0) will differ from previously published "virial" 
results In a second step, we will derive the next correc- 

tion to SC-0 in the str ong -coupling expansion, following the 
Wigncr picture of 0, [l6[ , since it has been shown that the 
virial approach fails in this task [l^|. Section HVl will contain 
the essentials of the Monte Carlo method used for numeri- 
cal simulations. Since the evaporation phenomenon exhibits 
particularly pronounced finite-size effects, we will resort to 
exponentially large system sizes [with typically \og{D/ R) on 
the order of a few hundreds], adopting the efficient centrifu- 
gal sampling scheme devised in Ref. The comparison 
between the analytical predictions and simulation data will 
be provided in sections [V] and IVD Some emphasis will be 
on density profiles, but an order parameter for the conden- 
sation phenomenon and pair correlation properties will also 
be analyzed (section [V]). While most simulations have been 
performed with system sizes that are large enough and not 
plagued by finite-size effects, the consequences of decreasing 
the system size will be addressed in section IVI1 Finally, sec- 

The technical content of 



tion IVIII contains our conclusions, 
the presentation has been lightened by delegating details to 
appendices. 



II. MEAN FIELD 

We recall in this section some known results from the 
mean-field (Poisson-Boltzmann) theory [12]. The non- linear 
Poisson-Boltzmann equation admits an analytic solution of- 
ten attributed to Katchalski et. al. \\SM1V{ . but which 



seems to date back to Liouville |22J. This solution brings 
to the fore the importance of a large lateral extension pa- 
rameter [A = \og(D/R) with R and D the charged cylin- 
der radius and the outer cylinder cell radius respectively, 
see Fig. [3]. Poisson-Boltzmann equation reads, outside the 
charged cylinder (r > R) 



V|u(r) 



(1) 



2 



where u(r) = — /3ge$(r) is the dimensionless potential, and 
r the radial distance. Tilde distances are made dimension- 
less with the Gouy-Chapman length p = R/£ (r = r/p, 
R = R/p = £), and kjj is a constant that has no sig- 
nificance before a reference value (a gauge) is chosen for 
the potential. Eq. ([1} is supplemented with the boundary 
conditions d?u(D) = (global neutrality of the cell) and 
R(dru(R)) = — 2£ (from Gauss' theorem, normal component 
of the electric field proportional to the surface charge of the 
cylinder) . 
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which determines the behavior of p. From 
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we get 



-a(A + 1) 



(7) 
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The analytic solution depends on the Fuoss critical param- 
eter £ c = A/(l + A). Here, 
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where a is given by the transcendental equations 



if € <£ c 
(2) 



1— a coth 
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The corresponding dimensionless density p 
reads 



^) = ~2 X 



sinh 2 (a log ~ + coth 1 ^ 



sin 



a log i + cot^ 1 ^ 



(3) 

p/(27rZ B <7 2 ) 

if e < e c 
if e > e c 

(4) 



with a normalization condition 



p(r)rdr = £. 



(5) 



A related quantity of interest is the counter-ion integrated 
charge in a cylinder of varying radius, which is, using Gauss' 
law, 



Mr) = , ru '(r) 
A 2£ 



r £ — 1 

1 + a cot ( a In h cot -1 

it a 



(0) 



Electro-neutrality imposes that X(D) = A while A(i?) = 0. 

The choice of units (2ttIb<j 2 ) to measure the density is 
of course not essential, but proves convenient in that it will 
make contact densities at r = R of order one. In addition, 
the contact theorem [23j imposes that in the limiting case 
of an isolated charged plate, the contact density is strictly 
fixed to unity: p = 1. For a given surface charge <r, the 
planar limit is obtained taking R — > 00, other parameters 
being kept constant. It thus corresponds to £ — > 00 but S 
fixed, a thick cylinder case indeed (£ 3> S 1 / 2 ). However, as 
far as mean-field is concerned, the difference between thick 
and thin cylinders is immaterial, so that we should soon 
check that p(R) —> 1 when £ — > 00 after having sent the 
boundary to infinity (A — > 00). The SC profiles to follow in 
section lHll are not endowed with the same property, since the 
functional forms of p strongly differ in the thick and needle 
configurations. 

For the most part, the interesting regime is that of £ > £ c . 
In the large A = log(D/i?) limit, we then have cot(aA) ss 



p(r) = 



1 R 



A + l 



. (9) 



which holds for £ > 1. Hence, the density in the A — > 00 
limit 0] 



p{r) 
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1 + ^-1) log 
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(10) 



It appears that for £ < 1, p = at all distances, which 
signals complete evaporation of counter-ions. In other words, 
a cylinder is only able to bind ions if the value of the Manning 
parameter is higher than unity. For £ > 1, we further have 



p(r)rdr 



-1, 



(11) 



to be compared to ([5]). The evaporated fraction of counter- 
ions is therefore l/£. For large although not infinite values of 
A, the phenomenon remains, although of course, a distance 
criterion is required to differentiate condensed from evap- 
orated ions, since normalization ((SJ) always holds. To this 
end, a convenient inflection point criterion has often been 
used 

SIM HI (see also e.g. mi for a related discussion with 
added salt): the integrated charge A(r) plotted as a function 
of log r shows an inflection point precisely where 



log-|=A 



which corresponds to A(r c )/A 
fective integrated charge of £ 
r c = R at £ = 1 and for £ > 1 



oA 



(12) 



1 — l/£ and renders an ef- 
1. It follows directly that 



r c A + l 

log — ss 

6 R 2 



1 - 



1 



A(£-l). 

The distance r c is often called the Manning radius. 



(13) 



For values of £ > 1 ions will be condensed closed to the 
surface of the cylinder thus creating a cloud of charge with 
cylindrical symmetry that will screen the potential that other 
ions farther away perceive. Then, for a sufficiently large dis- 
tance -beyond r c — the integrated charge of the ion cloud 
and the cylinder is effectively £ = 1 (Manning Condensa- 
tion) and the ions in the outer region will interact with an 
effective cylinder with charge equivalent to £ = 1, which lies 
at the borderline of condensation. From the previous analy- 
sis, the fraction of ions condensed within r c is /m (Manning 
condensed fraction of ions), with 



hi — 



£-1 



(14) 
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and the condensed ions will form a cloud of charge £ — 1. 
From (fT5|) it appears that the Manning radius r c , for A very 
large, is close to A/2 in logarithmic scale, and the normal- 
ization condition reads 



p(r)rdr = / M £ = £- 1, 



(15) 



to be compared to (fTTj) valid for infinite dilution. 

On the other hand, to quantify the extension of the electric 
double-layer, it is also instructive to compute the distance 
at which the integrated charge is half the condensed charge 
(i.e. A(r)/A = (£ - l)/2£). From eq. ©, 



log 
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(16) 



This teaches us that the relevant length scale for the exten- 
sion of the ionic profile is R, a much smaller scale than the 
Manning radius. The above expression is compatible with 
the known fact that in the planar limit, the extension is 
given by the Gouy length p. Indeed, when R — > oo so that 
£ — > oo, we get from (fT6|) that x^^ — R oc R/£ = p. 

For future comparison with simulation data, we also pre- 
cise the behaviour in the vicinity of the outer cylinder. From 
eq. 0: 



1-log- 



(17) 



Conversely, in the vicinity of the charged cylinder when the 
Manning parameter approaches the critical value - i.e. £ — > 
£ c , we see from the transcendental eq. §5§ that for £ = £+, 
a — > thus yielding 



1 



A + 1 - log ■ 



(18) 



which gives a non- vanishing value for the density at contact 
of A -2 for the density. A similar approach to the one per- 
formed for £ > £ c renders the same result for £ = £~. Note 
that the value at contact is strongly dependent on the log of 
the box size. 

In the following analysis, the results displayed do not de- 
pend on the value chosen for the valency q of counter-ions, 
which will therefore not be precised, since it only matters 
through £ oc q and Sag 3 . 



normalized for £ < 1 -we consider here the infinite dilution 
limit-, yet another illustration of complete evaporation: the 
omitted prefactor is vanishing for £ < 1, so that po{r) = 0. 
For £ > 1, the previous profile is normalizable though, and 
wc have 
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(19) 



where wc have assumed a fraction of condensed ions /, i.e. 

/ p(r)rdr = ft. (20) 

We have seen in section [II] that / = Jm = 1 — l/£ within 
mean-field, a result that nevertheless holds beyond mean- 
field @. Taking / = f M = 1 - we write 



Po(r) 



2(g - l) 2 (R 

e 



SR^ 2 
r 



(21) 



This is our leading prediction, denoted SC-0 in the remain- 
der, which turns out to differ from the result derived in [H, [s| 
where the same form as (jTSJ) was considered, but with the 
choice / = 1 that turns out to be incompatible with evapora- 
tion of a non vanishing fraction of ions. We will refer to the 
choice / = 1 as the SC-0 / = 1 form. It can be noted here 
that the typical distance into which the ions are confined is 
given by R, and does not depend on 2. More specifically, 
if we compute the distance corresponding to confinement of 
50% of the ions, we get xfj^ = R2^~ 2 ^ , which exhibits a 
similar form as its mean- field counterpart (|16p . The coupling 
parameter 3 is indeed absent, but the expression is compat- 
ible with the ground state requirement that Xy 2 — > R when 
temperature vanishes, because then £ — > oo. For £ of order 
1 but larger than 1 to avoid complete evaporation, the rele- 
vant confinement scale is R. When £ becomes large, we have 
noted that x\/ 2 — R^ 0, and more precisely Xi/ 2 — R— R/£,, 
which, again, is the Gouy length fi oc (qlga), setting the con- 
finement range in the planar case. 

We note that the contact density following from (|21[) reads 
Po(R) = 2(£ — l) 2 /£ 2 , that is exactly twice the mean-field 
contact density found in (jTUJ) . It also appears that it is not 
possible to recover the planar limit with its p(R) constrained 
to unity by the contact theorem, since expression (|2ip only 
holds for small values of £ 2 /S -this is the needle constraint-, 
while the planar limit is met for £ — > oo at fixed S. 



III. STRONG COUPLING 
A. Leading order behaviour (SC-O) 

In the strong coupling large 3 limit, at fixed £, the ratio 
a' /R oc S/£ 2 becomes large: this is another way to define the 
needle limit. Hence, the typical distance between particles 
becomes large compared to the radial distance they explore, 
and to leading order the same single particle picture as in 
the planar geometry [13|, [ljj does hold @ • The ion-ion inter- 
actions become sub-dominant compared to the rod-ion term, 
and the counter-ion profile is thus given by the exponential 
of the bare cylinder logarithmic potential. This means that 
Po{r) oc r -2 ^ The subscript refers to the dominant 

order in a large 2 expansion. This functional form cannot be 



B. Correction to leading order (SC-1) 

Before comparing our SC-0 prediction to numerical data, 
we adapt the method used in Refs. [ill 13 to compute the 
next term in the strong coupling expansion. Such a pro- 
cedure yields fundamentally different results than the virial 
approach of Refs. [l7|, [13, since correction terms appear 
dressed with a different power of S. These corrections have 
nevertheless not been worked out at virial SC level in the 
present cylindrical geometry. 

The starting point is to determine the ground state of 
the system, and to further consider the relevant excita- 
tions, those which contribute to the correction to SC-0. In 
other words, we should identify the excitations from parti- 
cles ground state that have the smallest energy cost, and 
we therefore first expand the inter-particle potential assum- 
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Figure 5: A R — >■ approximation of the ground state for the 
cylinder system in the needle limit. The length a' should 
account for ionic evaporation, so that it is defined as 
a' = q/(\f), with / = 1 - Hence, 
R/a' = C 2 //H = £(£ - 1)/S. 



ing that the particle displacement X (X: = x + zfc & 
X 2 = a; 2 + z 2 , x being a vector in the plane perpendic- 
ular to the z axis) from its lattice position is small, i.e. 
|X| -c a'. We will assume that ground state positions are 
given by R; = a'lk, for I 6 Z, which corresponds to parti- 
cles localized onto the z axis, see Fig. [5l This is consistent 
with the needle limit where R/a' — ¥ for S — > oo. A more 
correct ground state is sketched in Fig. [3J but for computing 
the desired correction to po, it is sufficient to distort the true 
ground state into the simpler form shown in Fig. [5j the same 
leading correction ensues. 



The energy cost for a given configuration of ions reads (see 
appendix [Al for more details), 
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(22) 



where £ is the Riemann Zeta Function (£(3) ~ 1.202). In all 
formulas and provided that the system size is big enough, 
we consider / = 1 — l/£. We however leave the fraction 
/ apparent, for it becomes a non trivial function of £, cou- 
pling parameter, and confinement when finite-size effects do 
matter (see section IVI[) . 



condition feq. Hip , we arrive at 



1 + C(3)4/ 3 



£-1 



r) ~ e-2 



(23) 



The form of Eq. ([22j) calls for some comments, before its 
use in the Wigner strong-coupling machinery fl3l fl(| . When 
considering separately the displacements of particles along 
the rod (variables z), or perpendicularly (variables x), a dual 
localization phenomenon appears. The coupling parameter 
S governs the localization of ions onto the rod, which was 
already clear from the relation R/a' oc £ 2 //S < £ 2 /S. It is 
not coupled to the z degree of freedom, so that even at very 
large 3, the essentially one dimensional system of condensed 
counter-ions may be fluid. It is then the Manning parameter 
£ that governs crystallization along the rod direction, and 
is thus the parameter coupled to the z degree of freedom in 
([2"2"]l . Another way to recover this conclusion is to compute 
the coupling (plasma) parameter corresponding to ions on a 
line, with inter-particle distance a': we get q 2 lB/a' , which 
is equal to £ — 1. We therefore expect a one dimensional 
transition for large values of £, at large S such that the needle 
scenario holds (S ^> £ 2 ). In all this discussion, it is implicitly 
understood that the ions are typically confined, radially wise, 
in a sheath of extension R around the cylinder. This is indeed 
the case, in the SC regime as well as within mean-field, see 
the discussion in section MI Al 



Considering the Boltzmann weight constructed from the 
energy (|2"2"1) . we fix one tagged particle at a given position 
Xo, and integrate over the remaining particles, in the spirit of 
the procedure worked out in Refs. [13|,[l6|: /o(xo) = C(<5(x — 
xo)), with C a normalization constant. After some algebra 
detailed in appendix ^ and under the proper normalization 



subsequently referred to as SC-1. We will discuss later the 
limit of validity of the above expansion. It should be empha- 
sized that the profile (|23p is an expansion in r, which does 
not hold up to oo. This poses a problem for normalization, 
since the neglected higher order terms become prevalent for 
£ < 2, and explain why (|23p can only be normalized for 
£ > 2. This leads to the conclusion that while the zeroth or- 
der term po may give a reasonable profile for small values of 
£ (to be precised in section |V| , the correction pi is deficient 
for £ < 2. 



C. A single particle variant 



There is an alternate semi-numerical treatment to the 
strong-coupling problem. Within the range of high values 
of £ and 2, such that £ 2 //S <C 1, we fix the ions at their 
ground state positions, and compute the energy cost if one 
particle (and only one) is shifted perpendicularly from this 
structure, consistent with the single particle picture in the 
strong coupling regime. The density profile is then 



p(r) oc ( — 
r 



2« 



exp 



j>0 



(24) 
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up to a normalization constant. The series involved in the 
calculation has no known closed form. For large distances, 



j>0 



1 



\/j 2 + x 2 



— lOE 



const. 



(25) 



meaning that the large distance behaviour is p(r) oc l/r 2 , as 
is the case within mean-field. The profile (|2~4")) is, therefore, 
not normalizable when A — > oo. There is however a large 
range of upper cutoff distances where the resulting normal- 
ized expression ((24)) is invariant close to the charged rod, so 
that the normalization problem can be in practice easily cir- 
cumvented. Eq. (|24p can be viewed as an improved version 
of SC-0, and will be referred to as SC-O*. In particular, it 
reproduces the r~ 2 ^ behaviour in the vicinity of the charged 
rod. 



IV. MONTE-CARLO SIMULATIONS 

For the numerical computations, we will adapt Monte- 
Carlo sampling to the cell geometry, and impose periodic 
boundary conditions along the main axis of the cylinder (z), 
taking due account of the long range of Coulomb potential. 
In the xy plane, the counter-ions are confined between the 
cylinder (R) and the outer shell (D). There has been exten- 
sive reports on numerical implementations for the calcula- 
tion of the potential for periodic boundary conditions under 
a number of geometries I2ja - l30l | . For the ID periodic case, 
the Lekner-Sperb sums [28[ have been the standard method 
to account for all the electrostatic contributions. The evalu- 
ation of such sums is numerically expensive, and requires in 
particular the calculation of an important number of terms 
for short distances. Here, we introduce a novel analytic for- 
mulation deduced from the Poisson-Jacobi transformation, 
as proposed in [3lj . It is simple to implement and is free 
of divergences in all ranges of interparticle distances. The 
resulting energy is based on the Ewald separation of the po- 
tential in a term that converges quickly in real space and 
another that converges quickly in Fourier space. Details arc 
presented in appendices [B] and [C] 

The potential energy of the system can be expanded as 



U = U R + U F + U c + U s . 



(26) 



Each of the terms are written in terms of two conveniently 
defined variables, p^ the distance between the particles' po- 
sitions projected to the plane perpendicular to the z axis and 
here L z is the length of the box along the the 
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periodic direction which is naturally chosen as the cylinder 
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where qi is the charge of particle i, a > is a parameter 
chosen for convergence, 7 the Eulcr-Mascheroni constant, 
Erfc(x) the complementary error function and Ei(x) is the 
exponential integral as defined in Appendix [C] 
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where k = 2im/L z for n € Z and Ko(x,y) is the incomplete 
Bessel function. Further references for the evaluation of this 
function can be found in [321. [33j . 
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with A the linear charge density of the cylinder (a = 27Tf)- 
Notice that the cylinder is located at the origin of coordi- 
nates. 
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(30) 



The previous expression for the energy is free of divergences 
for any value of = and since the particles are bounded 
to the cell, we have < \og{pi/R) < A. All simulations were 
taken with a number of particles that ranged between 300 
and 1000, and ran typically over 10 7 steps. 

In the problem under scrutiny, finite-size effects are impor- 
tant and logarithmic in D/R [9j, which requires very large 
system sizes and consequently precludes standard sampling 
methods. To circumvent this difficulty, we use the centrifugal 
sampling technique 0] which consists in writing the partition 
function with more convenient log variables, Y = \og(r/R). 
Then, the partition function 



Z = Co I dv N dz N exp (-f3U(r N ,z N )) (31) 



v N 



transforms to 



Z' = C, 



I J d6dY N dz N exp (^-(3U(r N , z N ) + 2 ^ Y^j , 

(32) 



which redefines the energy into 



U<(Y N ,8 N ,z N )=U(r N ,z N )-±Y, Y >. (33) 



The calculation of the inter-particle potential energy requires 
to know the Cartesian coordinates of the particles; hence, the 
transformation of coordinates has to be performed each time 
a particle moves. For the latter choice of variables we can 
choose a Monte-Carlo step size of (AY, A9, Az) such that for 
short distances (7 Ri 0), a displacement of Ar is at most of 
a Gouy length. Equivalently, 



AF = los 



R + p 



R 



(34) 



For the sake of efficient equilibration, it proves useful 
to distinguish between two kinds of ions in our system. 
Bounded ions arc constrained to a shell measured in R 
units as discussed in sections |TT] and MI Al Their unbounded 
counter-parts cover a region beyond the Manning radius (r c ) 
where the average inter-particle spacing is much larger than 
the electrostatic correlation length Zf,, hence forming a weakly 
coupled gas. From mean field (eq. (H3J)), we expect the lo- 
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cation of the Manning radius near A/2 in log units. This 
question will be further explored when presenting our re- 
sults. As a part of the equilibration process, particles have 
to be exchanged between the two populations. To this end, 
we propose to move one particle from Y to Y' = A — Y. 
Performing twice such a move will return a particle to its 
original position. Therefore, detailed balance is preserved by 
choosing a fixed probability p (p ~ 10~ 4 ) to select this type 
of particle exchange move over a regular one. The attempt 
is then accepted employing the usual Metropolis criterion. 
Such an approach guarantees proper equilibration indepen- 
dently of any chosen initial condition for the counter-ions. 
This is particular important for the study of finite size sys- 
tems at large coupling, as algorithms with only standard 
moves (including the original centrifugal sampling method) 
will not sample correctly the configuration space. 



V. PROFILES, CORRELATIONS AND ORDER 
PARAMETER 

We present in this section the bulk of our results. A large 
system has been simulated, in order to get rid of finite-size 
effects, that shall be studied separately but in a more cur- 
sory fashion in section [VTl We start by validating our sim- 
ulation procedure against known results. To this end, two 
features can be used. We can first check that for small S, 
the mean-field (MF) expressions are recovered, and also that 
at arbitrary 5, the condensation scenario coincides with the 
MF one §. 



A. Counter-ion condensation / evaporation 
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Figure 6: Monte Carlo measured cumulative density of particles 
for A = 300 as a function of the logarithmic distance log(r/_R), 
for H = 10 _1 , 10 2 , 10 5 and £ = 1, 2, 3, 4. The crosses drawn close 
to (A/2, / = 1 — l/£) are the expected locations of the inflection 
point, at the Manning radius r c given by (|13[) . The mean-field 
prediction (|6j) superimposes to the H = 0.1 results. 

We begin by the condensation phenomenon. For £ < 1, all 
ions dilute away when A = log(D/i?) increases. For £ > 1, 
a fraction / of ions remain condensed in the vicinity of the 
charged rod. This can be seen in Fig. [SJ which shows the 
integrated line charge A(r) in a cylinder of varying radius 
r, as introduced in the mean-field section [TTJ By definition, 
A(i?) = while electro-neutrality requires A(Z?) = A. The 
ionic atmosphere is more bound to the rod as the coupling 
parameter S increases, and the profile then strongly departs 



from MF. Further from the rod, ionic correlations decrease 
as a consequence of the lower ionic density, to such an extent 
that the tail of the ionic profile is described by mean-field. 
This is why in the right hand side of the figure, the different 
S-curves collapse, and coincide with MF form. The inflec- 
tion point property, that is clearly visible, hence takes place 
at a point that is EE independent, for which mean- field re- 
sults apply @. This provides the rationale for the two- fluid 
picture (bound ions before the inflection point, unbound be- 
yond) that is often used for polyelectrolyte, but that is quite 
specific to the salt- free case p. Note also that to observe 
the MF condensed fraction / = /m = 1 — 1/6 as in Fig. [51 
exponentially large box sizes are required. We come back 
to this point in section PVTl Our results reproduce previously 
reported data [jl|9|. For instance, extracting the inflection- 
point location from plots such as Fig. [51 we always obtain a 
condensed fraction that is extremely close to /m = 1 — 1/6, 
see Fig. [71 Figures [HI and [7J justify the normalization choice 
made in section HIT1 that led to Eqs. (j2l"j) and ((23)) . Indeed, 
the strong-coupling profiles are meant to describe the ionic 
atmosphere in the vicinity of the charged rod, an atmosphere 
that is deprived, over an exponentially large distance range, 
from the ions that lie in the vicinity of the confining border. 
This results in the plateau of Fig. [51 at a value that does 
not correspond to full neutrality. Note that the term "vicin- 
ity" here should be taken in the broad sense, since it can 
be seen in Fig. \6\ that A(r) changes to reach full neutrality 
A(D)/A = 1 in the range where 260 < log(r/R) < 300, so 
that r changes by a factor e 40 ~ 10 17 . 
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Figure 7: Condensed fraction of ions as a function of the 
Manning parameter £ for values of S in mean field and in the 
strong coupling regime (A = 300). 



B. Density profiles 

Having validated our normalization procedure from inte- 
grated profiles, we turn with Fig. [8] to a more precise analysis 
of the profiles themselves close to the charged cylinder. As 
expected, the small EE results coincide with their mean-field 
form, see the EE = 0.1 curve. On the other hand, when EE 
is large enough and exceeds 10 2 , all profiles collapse onto 
the SC-0 prediction, Eq. (|2"Tj) . As anticipated, the mean- 
field contact value [p(R) = 9/16 ~ 0.56 for £ = 4] is half 
its strong-coupling counterpart. A more thorough analysis 

in 



of the contact density will be presented in section IV F 
conjunction with the study of the corrections to SC-0 and 
the test of the SC-1 formulation. 
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Figure 8: Radial density profile for different values of H and 
fixed £ = 4 for A = \og(D/R) = 300. The symbols show the 
Monte Carlo results. The lines are for the the mean field (MF) 
result and our SC-0 (strong coupling to leading order) analytic 
solution (|21[) . The SC-0/ = 1 prescription of Ref. [H,iJ is shown 
by the upper dotted curve. The reduced density is defined as 
p = P /(2txI b g 2 ). 



comes as a surprise that we can get semi-quantitative agree- 
ment for such low values as £ = F2, and good agreement for 
£ = 1.4. 
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Figure 10: Same as Fig. [5J where the SC-0/ = 1 prediction of 
Ref. 0] is shown by the dotted curve in the inset 



The algebraic form of SC-0 is better appreciated in the log- 
arithmic plot of Fig. [SI which shows that the 3 = 10 2 points 
tend to depart slightly from SC-0 for r > 1.6R, while those 
for 5 = 10 3 are in excellent agreement with the prediction 
([21]) . As explained above, Eq. ([21"j) is a double expansion, 
first in large 3, and second in small distances. It therefore 
does not hold up to arbitrary large r, and we will comment 
further the large distance behavior in section IV CI 
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Figure 9: Log- log plot of the radial density profile for different 
values of 3 with fixed £ = 4. The symbols display the Monte 
Carlo results. The dashed line indicates the mean field (MF) 
prediction coinciding with the S = 10" 1 Monte Carlo data. The 
continuous curve is for the SC-0 strong-coupling to leading order 
analytic solution (|2ip which superimposes to the 3 = 10 3 data 
and partially with the 3 = 10 2 points, and finally the dotted 
curve is for the SC-0/ = 1 formula of [5|, 

Consistently with the integrated charge plateau of Fig. [6l 
the / = 1 normalization of Ref. 0, Q does not fit with the 
numerical data, see the upper dotted curve in Fig. [5J This 
phenomenon is all the more pronounced as the Manning pa- 
rameter is low, see Figs. fTOl andfTTl On the other hand, the 
SC-0 Eq. ([2"T]) for / = 1 — l/£ provides a reasonable profile 
when 3 is large enough, and £ not too close to unity: the 
agreement in Fig. [TT]is correct, and better than in Fig. 
This was expected, since our SC approach is an expansion in 
the vicinity of the ground state of the system, and therefore 
better when, in addition to 5, £ is large enough. It even 
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Figure 11: Same as Fig. [TDJ for £ = 1.4. 



C. Universal crossover to mean- field at large distances 

The decay of the ionic profile with distance r has led 
various authors to surmise that mean-field should hold at 
large enough distances 0, H3, . Indeed, one may define 
a local coupling parameter S(r) from the ratio of Bjerrum 
length to the typical distance d between counter-ions at a 
given distance r. With p(r)d 3 of order unity, we expect 
3(z) oc p(z) 1 / 3 , which should bring the profile into the mean- 
field region S(z) <C 1 when z — > oo. We emphasize that 
explicit checks of this expectation are in general difficult to 
perform, due to the fact that the distances one should be able 
to probe can be very large. Our present study is nevertheless 
particularly well suited for investigating such an effect, due 
to the large systems considered. 

We start by considering infinite system sizes A — > oo. We 
note here that if mean-field holds at large r, then, the profiles 
should not only become independent of 3, but also on £ 
(attention should be paid here to the difference between the 
normalized profile p and the original one p cx Ib<J 2 J)). This 
is a property of the relation p0[) which gives asymptotically 
that p(r) oc £ -2 (i?/r) 2 / log 2 (r/i?). Going back to the initial 
profile, we get p{r) oc / i3 1 r -2 / log 2 (r/i?), for all values of £. 
This universality is illustrated in Fig. 1121 where the curves 
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for different £ and for large couplings asymptotically coincide 
with the MF expression. The figure also highlights the fact 
that the distances needed to evidence the MF form are quite 
large, r/R > e 10 ~ 2 x 10 4 . The figure corresponds to A = 
\og(D/R) = 300, so that for the range of distances displayed, 
the behavior is very close to its A — ¥ oo limit. A single mean- 
field curve appears in Fig. [T^l since for the range of distances 
shown, the MF solution for the different £ values differ only 
for very small r. 
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Figure 12: Plot of £, 2 p/q 2 = 2ttIbR 2 P as a function of radial 
distance, for A = 300. As in other figures, the symbols are for 
the Monte Carlo data. The mean-field prediction, which holds 
for the three values of £, has been added (line). 

Figure [T2] revealed that the mean-field tail is visible for 
small densities / large distances only. On closer inspection, it 
appears that the departure from the SC behavior which holds 
at small r is quite sharp, as shown in Fig. 1131 the algebraic 
profile in r~ 2 ^ holds up to r/R ~ 1.4. It abruptly evolves 
into another form for larger distances, a form that is still far 
from the MF expression (upper dotted curve) but closer to 
the critical £ = 1 mean-field curve (which of course is fully 
compatible with the MC results at the low 5 = 0.1). The 
latter remark provides an approximate means to compute 
the crossover point r c where SC-0 ceases to hold: we simply 
equate the SC-0 (eq. (j2"Tj0 and MF-£ = £ c density value at 
contact (eq. (|18p) forms to get 



logx c = - log \/2/jvfA 



(35) 



where x c = r c /R. We learn here that the dominant form of 
the crossover point behaves as 



As 



(36) 



The dependence on the Manning parameter agrees qualita- 
tively with Fig. [T^] where x c decreases upon increasing £. 
One interesting trait of finite A on the distribution of ions 
is that at the critical Manning transition parameter (£ = £ c ) 
the value of the contact density at the surface of the cylin- 
der is not zero but A -2 , which can be a small quantity, thus, 
providing the order of magnitude of the density that must 
be reached before a mean-field like behavior can emerge. 

A complementary means to illustrate the universality of 
MF behaviour at large distances is provided in Fig. [T31 which 
is explicitly governed by finite size effects. We know from 
section [TT] that the density in the vicinity of the confining 
cylinder at r = D is given by Eq. (fT7|) . The data collapse 
displayed in Fig. Q3]is remarkable, and shows that the £ and 
5 independent mean-field physics is at work in the tail of the 
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Figure 13: The radial density for different values of S, £ and 
A = 300. The dashed lines indicate respectively the mean field 
(MF) profile for £ = 1 and the dotted curve for £ = 25. The 
strong coupling to leading order analytic solution (SC-0) is 
presented in the solid curve. The right arrow indicates the 

crossover location (x c « 1.38) and the up arrow indicates the 
location of a Gouy length from the surface of the cylinder. 
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Figure 14: Profile near the edge of the box, for £ = 2 or 3, 
A = 300, and values of S in all regimes. The solid curve 
corresponds to the mean- field Eq. (11711 , 



D. Order parameter for the evaporation transition 

The mean inverse distance is a parameter that can con- 
veniently be used to see the transition from the condensed 
to the de-condensed phase for which, a priori, we expect a 
critical change around £ = 1. We therefore define the order 

parameter Si as A^ 1 X)i=i where the overline refers to 
the Monte Carlo time average at equilibrium. In other words, 
we have 



(37) 



Discarding box size effects (assuming A —> oo) Si behaves 
as, 

★ Mean field (MF - eq. (JTSJ)) 



S[ 



{MF) 



£ - 1 - e^Ei 



(38) 
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* Strong coupling to leading order (SC-0 - eq. ([15])) 

Msc-o) = f 2(C - 1) = 2(g- l) 2 
7 £(2£-l) C 2 (2C-1)' 



(39) 



Figure [T5] shows the results from the Monte-Carlo simula- 
tions taken with values of £ and S in all ranges. Si vanishes 
for £ < 1 since all ions dilute away from the charged rod, 
while Si ^ when £ > 1. As before, the agreement with 
mean- field is excellent at H = 0.1, and equally good is the 
consistency with SC-0 for large S (10 3 , 10 4 and 10 5 ). The 
data at 5 = 10 2 appear quite close the the strong-coupling 
limit, but exhibit some discrepancy. These data further illus- 
trate the relevance of normalizing the profile with / = 1 — 1 /£ 
and not / = 1, see the upper dotted curve. 
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Figure 15: Order parameter Si (mean inverse distance from the 
cylinder) as a function of the Manning parameter f , for different 
couplings S and A = 300. The lines indicate the mean field 
[MF, Eq. ([38"|l] and strong coupling to leading order [SC-0, 
Eq. (|39p ] analytic solutions. The upper dotted curve 
corresponds to the analytical prediction, Eq. (68), in Ref. @, [||. 



E. Towards the ground state: crystallization scenario 
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Figure f 6: Normalized pair correlation function along the 
cylinder for £ = 3 and different values of S. The number of 
particles is N = 300 unless specified in the legend, and A = 300. 



Upon increasing £, crystallization occurs along the cylinder, 
as hinted in Fig. [T71 where we recover the length scale a' as 
the correct measure of inter-particle distances along z. 
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Figure 17: Pair correlation function along the z-axis for 
different values of £ and S = f0 4 , and A = 300. The number of 
particles is TV = 300 unless specified in the legend. 



After having investigated the properties of the density 
profiles perpendicular to the charged rod, we now address 
the question of the correlations along the rod (z direction), 
through the computation of the essentially one dimensional 
pair correlation function g z of bounded particles along the 
z-axis. The normalization of this object was done with re- 
spect to the number of ions close to the surface of the cylin- 
der, thereby considering in the calculation only those ions 
present between two concentric cylindrical shells of at r = R 
and r = R + 10/i. Normalization ensures that g{z) — > 1 at 
large z. 

The results for g z are presented Fig. [TO] and [TTJ They 
shed light on the crystallization phenomenon that takes place 
here, and on the asymmetric roles played by the two param- 
eters 2 and £. Fig. [TrJl shows that beyond a certain value, the 
correlations along the cylinder no longer depend on 5. This 
was expected, from the dual localization argument developed 
in section IIIIBI A large S confines the ions in the vicinity 
of the charged rod, but their interaction along z is governed 
by £. If the latter quantity is small, we face an effective 
one dimensional liquid, that is mildly modulated in Fig. If6l 
In passing, this figure illustrates that the number of parti- 
cle taken for the simulations (300) is actually sufficient for 
our purposes: identical results are obtained with N = fOOO. 



From the form of Eq. (|22p and in particular the harmonic 
energy term in z, we expect that the first peak of g z will 
present approximate Gaussian shape. Under this assump- 
tion, the width of the peak can be read directly in Eq. (|22| , 
to be 

5 z ^{U)~ h - (40) 

From Fig. If 71 we can extract S z performing a Gaussian fit 
of the first peak. The resulting width is shown in Fig. [T8l 
which gives credit to the naive estimation and show a very 
good agreement with the linear trend expected for S z as a 
function of 8 Z cx (C/)" 1 ^ 2 - 



F. Ion Profile Close To The Cylinder and correction 
to SC-0 

We so far focused on quantities that were accurately de- 
scribed by the SC-0 form at large 3. Our goal is now to test 
the validity of improvements over this leading form (SC-1 
or SC-0*). To this end, we plot in Figs. [H [20] and ED the 
quantity r 2 ^p(r), that yields a horizontal line at SC-0 level, 
which is a way to make deviations from SC-0 more appar- 
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Figure 18: Width of the first peak of the numerical data for g z Figure 21: Same as Figs. [H and [H] for S = 10 4 . The solid 
in fig. [T7] The simple argument giving ((40]) leads to expect a curves represent the SC-0* prediction of Eq. ([24]). 

straight line. 



5- 1.5 
1.4 

(M 

S- 1.2 
1.1 
1 

0.9 
0.8 



J , 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

i ■ ^=2.8 


1 1 1 1 1 1 1 , 1 1 1 1 1 1 1 , 1 1 1 1 1 1 1 1 1 1 1 1 , 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 l_ 


I ■ £=3.2 




I ■ £=3.6 




\ £=4.0 













1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5 

r/R 

Figure 19: Radial ionic density for different values of 5, 
= 10 2 and A = 300. The dotted lines represent the analytic 
SC-1 result, Eq. ([231- 
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Figure 20: Same as Fig. 1191 for the same confinement but a 



higher coupling parameter 



cnt. The numerical data thereby obtained show an increase 
with r, compatible with a parabola, which is also the trend 
predicted by the SC-1 expression, see Eq. (|2"5|) . However, ob- 
taining a quantitative agreement requires considering large 
values of £, for the agreement displayed in Fig. is quite 
poor. The situation is better in Figs. |2"0"1 and |2"T1 We there- 
fore come to the conclusion that the leading SC-0 behavior 
may hold for fairly low values of £ as discussed in section 
MI Al while upon close inspection, the refinement SC-1 re- 
quires £ to be large. This comes as no surprise since ionic 



correlation do not enter the SC-0 form, while they are at 
the root of the SC-1 expression, derived assuming that all 
ions lie in the vicinity of their ground state position. As we 
have seen in section IV El this requires typically £ > 40, and 
explains the poor agreement in Fig. 1191 while we have better 
consistency in Figs. l2Tjlandl2n In addition, we have reported 
in Fig. [21] the results of the alternative improvement SC-0* 
as given by Eq. (|2~4"|) . It seems that such a route improves 
upon SC-0, but also upon SC-1. However, some care is re- 
quired in interpreting the results: while SC-1 follows from an 
exact although pcrturbative statistical mechanics treatment, 
SC-0* remains at the single particle level, and is heuristic. 
Hence, SC-1 may be viewed as providing the next to leading 
contribution in the SC expansion of the ionic profile, which 
is not the case of SC-0*. 

Note that if £ is too much increased at fixed 5, the needle 
requirement £ -C 3 1 / 2 may be violated at some point. This 
is not the case though with the data displayed in Figs. [121 
fJOlandini] For instance, we have in the worst case £/S - 5 = 
0.4. It should also be noted here that for large values of the 
Manning parameter, (r/R) 2 ^ becomes quite large for r/R> 
1; hence, small fluctuations in the density profile induce large 
fluctuations in the graphed results. 

To conclude this section, we report in Fig. [22] the Monte 
Carlo measures for the density profiles at contact, together 
with the SC-0 and SC-1 predictions. While SC-0 cxpectedly 
gives the correct main trend of p(R), it is seen that the 3 
dependent fine structure is well captured by SC-1. From the 
contact theorem, we also know that at fixed 3, increasing 
further £ ultimately leads to p(R) = 1. Such a trend is not 
visible in Fig. 1221 since the parameter range pertains to the 
needle limit, with constraint £ <C 3 1 / 2 . Additionally, if the 
limit 5 — > oo is taken first at arbitrary £, the SC-0 form 
becomes exact and we have p(R) = 2(£ — l) 2 /£ 2 , twice the 
mean-field expression as already noticed. Increasing next £, 
we finally get the large £ result p(R) — > 2, exactly twice 
the planar result. It is this trend that is illustrated in Fig. 
I2"2l To summarize, the limits of large £ and large S do not 
commute, and we can write 

lim Jim p(R) = 2 



lim lim p(R) 



1. 



The latter equality may be written more generally as 
lim^_ i . 00 p(R) = 1, for all 3. 
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Figure 22: Contact density versus Manning parameter. As 
above, the symbols are for the Monte Carlo data and the lines 
for the analytical predictions. The dashed line displays the SC-0 
form, which does not depend on S. On the other hand, the SC-1 
result following from Eq. (|23[) is 3 dependent, and there are 
therefore three different branches (continuous curves) showing 
the corresponding expectation for the three values of S studied. 



VI. FINITE SIZE EFFECTS 




log (r/R) 

Figure 23: Cumulative density of particles for A = 10. The 
inflection point is materialized by an arrow. The value of the 
condensed fraction, as given by the inflection point criterion, is 
indicated for each curve. 



In the previous section, we reported results for large sys- 
tem sizes, in order to have a universal fraction of condensed 
ions, and a clear cut distinction between those ions that par- 
ticipate in the screening of the charged rod, and the de- 
condensed one that lie in the vicinity of the confining border 
at r = D. These de-condensed ions, which exist for all finite 
values of £, are the precursors of the ions which dissolve away 
when the system size is increased to infinity. For smaller sys- 
tems though, finite-size effects are very pronounced. This is 
illustrated in Fig. [23] which differs significantly from its 
large A counterpart, Fig. Although A = 10 in the fig- 
ure is large enough to allow for discriminating condensed 
from de-condensed ions, it is seen that the inflection point, 
which still provides a convenient cut-off for the partition- 
ing, severely depends on the coupling parameter 3. In the 
mean-field regime, the corresponding condensed fraction is 
still given by / = /m = 1 — (hence 2/3 on the figure). 
Beyond mean-field, the quantity increases with S in a non- 
trivial fashion. 

A similar conclusion holds concerning the system size de- 



pendence, see Fig. [5H For A > 50 only do we get a con- 
densed fraction that is close to its infinite dilution expression, 
1 — 1/^ = 2/3 on the figure. A valid question is then to see 
if the strong coupling prediction SC-0, Eq. (fT9|) . holds for 
the profile, with proper normalization / following from the 
inflection point rule. We can conclude from Fig. [55] that 
this is indeed the case: the Monte Carlo data are in good 
agreement with our prediction, where the only a priori un- 
known is /, taken from Fig. [53] The profiles are sandwiched 
between the limiting forms having / = = 1 — l/£, which 
appears to be a lower bound for the condensed fraction, and 
/ = !• 
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Figure 24: Build-up and size dependence of the integrated 
charge plateau for different values of A. Here £ = 3 and S = 10 2 
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Figure 25: Density profile for £ = 3, A = 10, against the 
/-normalized strong-coupling prediction. The values of the 
condensed fraction /, used in the SC-0 form Eq. (|19[) . are those 
which are read in Fig. 1231 

The remaining and final task is to quantify the dependence 
of / on system size A and coupling parameter, given that 
/ ~ fu when A is big enough, see Fig. [7J To this end, it 
seems that one should distinguish the regimes of low S where 
one essentially finds the mean- field result Jm, see Fig. l26l 
from the more strongly correlated cases. Figure [5JJ shows 
that in the latter case and for fixed coupling S, / decreases 
with confinement in such a way that 



/-/. 



A I 



e 



(41) 



with e and 7 two dimensionlcss parameters reported in Table 
[H and obtained from regressions on all data sets collected. 
A first conclusion which can be drawn is that except for too 
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Figure 26: Condensed fraction of ions determined numerically 
using the inflection point criterion in the mean field and weakly 
coupled regime for £ = 3. 
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e 




7 




10 2 


3 


12.5 


2.2 


1.99 


0.08 


10 3 


3 


7.52 


0.06 


1.198 


0.003 


10 4 


3 


16.0 


0.8 


1.23 


0.02 


10 2 


4 


2.5 


0.2 


1.57 


0.05 


10 3 


4 


3.5 


0.2 


1.02 


0.02 


10 4 


4 


7.6 


0.5 


1.06 


0.02 


10 2 


5 


0.16 


0.13 


1.0 


0.2 


10 3 


5 


4.1 


0.1 


1.096 


0.007 


10 4 


5 


8.4 


0.2 


1.100 


0.005 



Table I: Numerical values for e and 7 from a regression 
performed considering — ^7 on multiple results. The 

quantities labeled a refer to the standard deviations. 



small A and £, we have 7 ~ 1. Quite expectedly, confining 
the system favors condensation. Second, the figure exhibits a 
departure from the 1/ A scaling on the left hand side, where 
A is small. In this region, the quantity plotted becomes A 
independent, and turns out to reach its maximum possible 
value, i.e. (£ — 1) _1 corresponding to / = 1. This is a hint 
that too small systems (say below some A;) do not allow for 
evaporation to set in. For A > A;, A(/ — /m) becomes A in- 
dependent, as also illustrated in Fig. [28l which furthermore 
shows that the 5 dependence is logarithmic. More precisely, 
we have 



/-/. 



M 



a(logH- 8) 



f. 



(42) 



with a and 6 dimensionless parameters reported in Table 
im All previously described effects are encoded within this 
ultimate expression. Note however that this result cannot 
hold for arbitrary large S, for / again has to remain smaller 
than unity. Hence, and similarly to small A results, full 
condensation with / = 1 is achieved at large couplings while 
holding a fixed size (see the £ = 5 and A = 30 last two 
points). 
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Figure 28: Same as Fig. 1271 to probe the 3 dependence in a 
log-linear plot. The dashed lines represent linear regressions 
that are summarized in Table HT1 Note that since / is bounded 
from above by 1 (complete condensation), the quantity plotted 
cannot grow without bounds which explains the departure from 
scaling at large 3. 
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<7S 


a 




30 


4 


4.42 


0.03 


1.37 


0.01 


50 


4 


4.45 


0.02 


1.276 


0.007 


60 


4 


4.42 


0.03 


1.228 


0.009 


300 


4 


3.9 


0.2 


1.01 


0.05 


30 


5 


4.50 


0.04 


1.25 


0.01 


50 


5 


4.55 


0.03 


1.22 


0.01 



Table II: Numerical values for a and 5 from a regression 
performed considering ^~^ K " — ^QSS£^°1 on multiple 
results. The quantities labeled a refer to the standard 
deviations. 
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Figure 27: Condensed fraction under strong-coupling, for 
£ = 3, 4, 5, in a log-log plot where errors are less than the tick 
size. The dashed lines represent linear regressions which are 
summarized in Table [I] 

Another interesting feature emerging from Fig. and 
Table |TT] is that the x-axis intercept (<5) is the same -within 



numerical accuracy- for the different sets (5 ~ 4.5). Con- 
sequently, for any given value of the coupling below e , the 
system will exhibit "ideal" evaporation (Manning evapora- 
tion or / = Jm) regardless of the size and, from the data, 
also regardless of the Manning parameter. This leads us to 
believe that it is a universal property in the evaporation of 
ions. Together with the fact that a takes values close to 
unity except for too small A, we summarize our finite-size 
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analysis with the expression 

f-fr ^ (logS- 4.5) 
/a/ A 

which holds provided 3 > e s ~ 90 and A > A;, while / ~ 
/a/ for E < e s . Full condensation determines A/ (at / = 1) 
in such a way that 

A,~(£-l)(logS-4.5). (44) 

It should be kept in mind that A = log(D/R) measures the 
log-size of the system. 



VII. CONCLUSIONS 

The behaviour of counter-ions at a planar interface is en- 
coded in the coupling parameter 3, defined from the uniform 
surface charge a as 3 oc l 2 B aq z . This scaling simply follows 
from the fact that the relevant Bjerrum length for q-valent 
ions is q 2 Ib, and that the natural measure of surface charge 
is a/q. Hence the dimensionlcss charge (g 2 /s) 2 cr/g. Equiv- 
alcntly, we can view 3 (or more precisely y/S) as the ratio 
of thermal energy over the typical Coulomb pair repulsion 
when all ions are condensed onto the plane, and where the 
inter-ion distance reads a oc dqjcr: 3 1 / 2 = (q 2 lB)/a. For 
large 3, the ions are confined in a region with extension given 
by the Gouy length \i oc (is erg) -1 . This applies to all cou- 
pling regimes, from mean-field at small 3 to strong-coupling 
[141 ] . When considering curvature, and addressing cylindrical 
macro-ions rather than planar, a new parameter enters the 
description, £ = qls^ = qlB^itRa. Of course, when R^> a, 
which also means £ 3 1 / 2 , the situation is very close to its 
planar counterpart. It changes significantly in the opposite 
case £ <C S 1 / 2 where curvature is strong, and that we called 
here the needle limit. There, the ions are confined in a region 
of extension R, given by the radius of the charged rod, that 
is S-indepcndent. In the needle regime, the radius R is much 
smaller than the typical distance a' between charges along 
the rod (R -C a -C o'). This provides the rationale for de- 
riving simple strong-coupling (SC) predictions: in the sense 
R <C a', the ions are far way from each other, and mostly 
respond to the log potential of the rod. A single particle 
picture holds, that can be seen as the ideal gas behaviour of 
non interacting particles in an external field (leaving aside 
the subtlety of counter-ion evaporation, that is a collective 
effect). Good agreement with Monte Carlo results can then 
be achieved, even at "small" values of £. This even if, strictly 
speaking, the ground state of the system is approached when 
both 3 and £ are large. The reason is that the single particle 
picture is already an acceptable approximation when £ is of 



order unity (or slightly beyond), provided 5 is large. We also 
add that again for large 5 where the ions remain close to the 
cylinder, the Manning parameter £ controls the essentially 
one-dimensional structure of the ionic system, from a repul- 
sive liquid at small £ to a crystal at large £. Indeed, although 
£ is initially defined as the dimensionlcss linear charge of the 
rod, it can be rewritten as q Is/a', and thus qualifies as a 
one dimensional plasma parameter, quantifying the strength 
of interactions for ions along the rod. 

An additional feature pertaining to cylindrical charged 
macromolecules is the counter-ion condensation-evaporation 
phenomenon, that plays a prominent role here. While its in- 
fluence on the mean-field behaviour is well known, we have 
shown that due account of its effect is essential for a good 
agreement between the SC theory and simulations. We have 
explicitly worked out the leading ionic profile under large 3 
(SC-0), together with the first correction (SC-1), that indeed 
improves upon SC-0, but requires quite large values of £ to be 
relevant (£ > 10). We also illustrated clearly that the ionic 
profile, even at large 5, crosses over to mean-field behaviour 
far from the plate. This behaviour, although expected, is in 
other settings extremely difficult to observe and could only 
be evidenced due to the (exponentially) large sizes used in 
the simulations. 

We have performed a finite-size analysis for the condensed 
fraction, which leads to several novel features. To this 
end, we proposed modifications to the previously introduced 
Monte Carlo sampling method, that significantly improved 
convergence rate. The signature of finite size effects is loga- 
rithmic both in the coupling parameter and in system size, 
provided 3 is above e s ~ 90. Full condensation can be 
achieved for both small A and large 3. 

Before concluding, we provide some parameter values for 
an important rod-like biopolymer. With double-stranded 
DNA, one has in water at room temperature (Is — 7 A), 
one has £ ~ 8, 5 ~ 22 with divalent ions q = 2, and £ ~ 16, 
5 ~ 180 with tetravalent ions. This latter case is not quite in 
the needle limit since ^/S 1 / 2 is of order 1, but approaching 
it. The expressions derived here, which are salt-free, would 
then provide a zeroth order limiting case. 

Interesting perspectives opened by this work include di- 
electric systems, the study of the effect of salt (added elec- 
trolyte), together with working out the two-dimensional pen- 
dant of our investigation, where ions interact with a log po- 
tential, which should lead to a large distance physics that is 
no longer of mean- field type (3^ . 
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Appendix A: Inter-Particle Potential Energy and Density Profile 



In this appendix we derive the density profile in the strong coupling and needle limit up to the first correction in the small 
"needle" parameter Rj o! = £ 2 //S. The following will be the notation used for normalization: x = xj [i and x = x/a', and 
bold letters for vectors accordingly. In this limit, the distribution of counter-ions presents small thermal fluctuations from 
the ground state. Then it is natural to write the position of a given particle as R; + X with R/ = a'lk and X = x + zk, 
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where |X|/|R| is of order R/a'. Then, 



1X.(2R + X) , 3(X.R)> Lop/fl?) 



|R + X| |R| 2 |R|3 2 |R| 5 

( 1 X • (2R + X) 3 (X ■ R) 2 
I ~2 jSp + 2 IRI5 



,i 2 _i^yi _i^yi_ i i (ex i B gH = j_?f 







i 1 




° i ii 




3 


1 i? e 2 


3 


1 r ~ hq 2 e 



-Rl 3 J 



where 

~a' 3 ~e\a'J R ~ e \a'J R ~ l B q 2 £ 3 \~ J J R ~ l B q 2 S 2 
Therefore, the ion-ion energy term can be written as Ylk<j V( a 'li — ^|) + with 

' 2S 2/ ^\ 2 IR.-RI 3 2 |R,-RP 

-- V - — —~ — ((R, - R) • x.,)((R 7 - R) • XA 
2^|R i -R ! |5 VV J j Jjl 

Using Rj = a'jfc, then Rj • X; = a'j(zi), 

mi f 3 v J ~ £ j 2 + 25 j v 1 i lv *j : *j + _ ^ v ? ' 

E 2-l 2 ^b'-Zl 3 2^ I7-ZI 3 2 ^Ij-ZI 3 

To evaluate the series we use the definition of the Riemann (^-function q) : = (q+j) n ano - : = 2^=1 



¥j u 1 ¥j u 1 

In order to fully analyze the crystallization effect we need to write the energy term separately for perpendicular and parallel 
displacements to the surface of the cylinder. The strength of the energy cost for displacements parallel to the surface have to 
be measured against the lattice constant a' (i.e. z = z/af), thus we have to write (36U in two separate terms re-normalized 
appropriately. Hence, 



m 4 f z\ -as,*? + 1 g fz| U if (f) 2 e Urn - * g ^ 



that yields the final result for the potential energy change 



^=i^\ + \ E fz| ) + U E { 2C(3) ^ 2 - * ■ g ITTF ^ (A1) 

Finally, the total energy, up to a constant term, including the cylinder potential term is as follows 



m =i £ log ^ J + fc/ 3 E | -c(3). 2 • iEy-y^j • -''E | 2 c(3)^ 2 - $ - £ ^ }■ ( A2 ) 

The particle density profile can be calculated as p(x) = C(5(x — xq)) (with C a normalizing constant satisfying eq. [TTj) . 
Notice that we have chosen the perpendicular displacements for the density profile. 



2 \ £3 



£ f3V^ J /-/ Q N~2 , 1 X Z -Xj I ^ JoA^s2 - \ ^ h 



p(±) =C (j[J dx fc j |^TJ y dzjj exp[-« , 



(A3) 
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where integral over the perpendicular displacements X& is over N — 1 particles. Expanding exp(—/35U) up to the first 
correction in R/a', we have 

with Co and D k normalization constants. The contributions from the z direction has been integrated out. We define x) 
a distribution as 

m,*): =£) (fr) > ( A4 ) 

where D = fci , V£ € 5ftA£ > 1 is a constant that normalizes the distribution. The domain of Q for x is ie e [R, oo) , 8 g [0, 27r]. 
To perform the calculations we will use the simplified notation 

(/) fc : = / dxi k n(R,x k )f, (A5) 



where (1)^ = 1. Using the latter notation, we notice that the probability distributions for different particles are uncorrelatcd 
((kj ■ Xfe) fc) j = ■ (xfe)fe, for k ^ j, and due to the cylindrical symmetry = 0. Furthermore (x|)fc = £ 2 (£ — l)/(£ — 2). 
Then, the density is 



p(x) =C (l + C('S)^f (V + (N- l)£ 2 f^) + o{{R/a'f)) ■ (A6) 

We proceed to evaluate Co such that / dv p(r ) = 27r/£, 



£2(i-£) ^ SW S 2J £-2 
where / is the ratio between the condensed integrated charge and the cylinder surface charge. The density profile is 

fep, =f«£f» (|)* (, + ^ _ izi) +o((fl/aT) ) , (AT) 

which yields a value at contact equal to 

Psc(R) =f 2 -^j^ (l - |^|i/ 3 + «((^/«') 2 )) • (A8) 



Appendix B: About the ID Ewald Summation 



Let us determine the Ewald summation for a one-dimensional periodic system. First, we write the potential as a summation 
of one particle interaction and its images at a point r (L„ = nLk for neZ) 



Y — 
E 
E 



r + L r , 



il e -|r+L„| 2 t 
I" Jo Vt 

I il p -|r+L„| 2 t 



(Bl) 



E 



1 ^ rft c -|r+L„| 2 t 
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with a a real positive parameter chosen wisely for convergence. The first integral is the well known Erfc(x)/x function. Then 
using the Poisson-Jacobi transformation for the second sum (exchanging sum and integral) 



E 

L„ 

E 
E 



Erfc(a|r + L„|) 1 f a dt 

|r + L n | ^ J Vi 

Erfc(a|r + L„|) 1 I"** dt 

|r + L„| 0rio Vi 



E< 



-|r+L„| 2 t 



E e " (z+? 



Lft 



Erfc(a|r + L n |) 



(B2) 



|r + LJ 



dt- 



Vt 



L 



E 



e 4t e 



ikz 



where k = 2 f' n for n £ Z. Simplifying and considering independently the term k = from the rest, 



^ Erfc(q|r + L n |) 1 \ - ikz 

2^ i r + L i + T M e 



fe^O 



* 



1 


' z 1 " 2 p-pV 




/ dt e 




7o * 



(B3) 



The second integral can be rewritten by substituting u = ^- = — *) 

Erfc(a|r + L„|) 1 



E 



Ir + L, 



fc#0 



e *° 



1 


2 2 , " 

r a p-pt 




/ di 6 


+ L 


Jo * 



(B4) 



where the second integral is by definition (|C6[) the incomplete Bessel function Ko(^-,a 2 /5 2 ), which can be represented with 
the uniformly convergent series [32l HH 



Ko(x,l/) = r(0,x + 2/)Jo(2y) + £] 



n=l 



r(n,x + y) 



(x + y) 7 



+ (-l) n (x + y) n T(-n,x + y) 



(B5) 



with r(n, x) and J n (x) the incomplete Gamma Function and the Bcsscl function respectively. The evaluation of the function 
for the y > x takes some time but we can use the following identity, 



Ko(x,y) = 2A' (2Vx^) - K (y,x), 
with Kq the modified Bcsscl function of order 0. Therefore, we can write the overall expression 





1 


' r a2 e~p 2t ~ 






/ dt 


+ L 


Jo t 



(B6) 



(B7) 



Taking e > 0, 



f/f- 



lim 
lim 



dt- 



dt- 



t 



dt- 



-P't 



t 



(B8) 



lim Ei(ep 2 ) -Ei(a 2 p 2 ) 

-7 - log(a 2 p 2 ) - Ei (a 2 p 2 ) + log a 2 - lim loge, 

e->0+ 



where Ei(x) is the exponential integral. Substituting 
Erfc(a|r + L„|) 1 



E 



Ir + L,, 



fe^o 



fc 2 

4a 2 ' 



iE^ (^aV] + T 



-7 — log(a 2 p 2 ) — Ei(a 2 p 2 ) + log a 2 — lim loge 



(B9) 



17 



The self energy term is found following the latter procedure excluding from the summation the L = term, 



^ ErfcHL^I 1 ^ / k 2 \ 1 



L„#0 1 1 fe^O v 

Erfc(a|L„|) , 1 ~ / k 2 \ 1 



log a — lim log e 

£->0 + 



1 eft 
^ Jo Vi 



E 

L„#0 



log a — lim log e 

£->0+ 



2a 



(BIO) 



Finally, we consider the contribution due to the cylinder within the cell and the particles and with its images. The following 
is the result of integrating with respect to z which leads the following two results, 



b cy i =V 1 (a 2 p 2 ) + 

= —7 — log(a 2 p 2 ) + loga 2 — lim loge 



-7 - log(a 2 p 2 )- Ei (crp 2 )+ log a - lim loge 

e->0+ 



(Bll) 



»cyl0 



= L E 



Erfc(a|L„|) 



E E: 



\4a 2 J 



log a — lim log e 



2aL 



L„^0 1 1 fc^O 

The overall energy [/ can be written as 



(B12) 



(B13) 



with A the linear charge density of the cylinder (a = )■ Notice that the cylinder is located at the origin of the coordinates. 
For neutral systems (J2 i qt + L\ = 0), 



U = U r + Ur + Uc + U s , 



(B14) 



with, 



JV-l N 



i=l J=i+1 



Erfc 



(« (4 + ( 



{pIj + (Zij + Ln ) 



if pij = 

2 s ) 3 ' £ 1— 7 — log (a 2 p 2 j) - Ei (a 2 p|j) i otherwise 



2 1 



JV-l N 



n E E E K ° (^2 > ) cos ( fc • 

u i=l j"=i+l fc>0 



(B15) 
(B16) 



(B17) 



47reen \ ^- i £n 

u \n>0 



Erfc(a-Ln) 1 ^ / k 2 \ a 
k>0 v 7 



(B18) 



Note that the term [log a 2 — lim e _j. + log e] in the overall energy cancels due to electro-neutrality. This final result is free of 
divergences and is absolutely convergent for all ranges of the variables pi , pij and 2y . 

I 



Appendix C: Special Functions 

The following arc definitions used throughout the paper. 

I The complementary error function Erfc(a;) 

poo 



I The incomplete gamma function T(a;, y), 

/>oo 

T(x,y) = / t x ~ x exp(-t)dt. 

J y 



(C2) 



Erfc(x) =-^= I e"* * 



' 7T 



2.r 



(CI) 



dt. 
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I The exponential integral Ei(x) t The incomplete Bessel function K„(x,y) 

00 (_i\k+l k r+00 e -xt-\ 

E 1 (x)^- 7 -logx- + ^ 1 L > X K v (x,y)= ~^ dt ■ ( C6 ) 

(C3) 



fe=i 



-di. 



Also, 

d 



. E x (i) = -— , (C4) 
dx x 



and the 2D Fourier transform is 

_ fc 2 

^[E 1 (aV)](k)=4 7 r 1 ~" — . ( C5 ) 
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